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The engineering tools of choice for the computation of practical engineering flows 
have begun to migrate from those based on the traditional Reynolds-averaged 
Navier-Stokes approach to methodologies capable, in theory if not in practice, of 
accurately predicting some instantaneous scales of motion in the flow. The migra- 
tion has largely been driven by both the success of Reynolds-averaged methods 
over a wide variety of flows as well as the inherent limitations of the method itself. 
Practitioners, emboldened by their ability to predict a wide-variety of statistically 
steady, equilibrium turbulent flows, have now turned their attention to flow control 
and non-equilibrium flows, that is, separation control. This review gives some cur- 
rent priorities in traditional Reynolds-averaged modeling research as well as some 
methodologies being applied to a new class of turbulent flow control problems. 
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1. Introduction 

In studying aerodynamic flows, new technologies in flow control, renewed inter- 
est in high-speed (sub-orbital) vehicles, drag and/or noise reduction and improved 
propulsion systems are now the driving research incentives. Although the usual 
Reynolds-averaged methodologies have well-known deficiencies, they do remain the 
front-line methodology for industrial prediction. Aerodynamic flows are character- 
ized by a wide parameter range associated with, for example, Mach number and 
Reynolds number and are turbulent or are transitioning to or from turbulent. Of 
course, the challenge from a practical engineering standpoint is to be able to pre- 
dict these flows so that design cycle times are reduced and designs optimized. It is 
clear that improving operating efficiencies are crucial in developing next generation 
aircraft and aerospace vehicles. For this reason, accurate and efficient prediction of 
turbulent flows has become an important topic. 

Even with the current level of computational power, the direct numerical simu- 
lation (DNS) of such complex turbulent flows is not feasible and will not be for the 
foreseeable future. The problem is simply the inability to resolve all the component 
scales within the turbulent flow. Even with computational grids exceeding 10 9 grid 
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points, simple fully developed channel flow simulations are still limited to Reynolds 
number (based on friction velocity) of O(10 3 ). Therefore, either all or part of those 
scales that cannot be directly computed need to modeled. The general task then 
is to develop means to partition the resolved and unresolved scales, and then to 
develop suitable models for them. Within the context of scale partitioning, two 
(related) approaches exist. 

The first approach is offered by the partitioning of the flow field into a mean 
and fluctuating part, an idea first proposed by Osborne Reynolds (1895) over a 
hundred years ago. This process, known as the Reynolds decomposition, assumes 
the instantaneous flow can be partitioned into a fluctuating part representing all the 
turbulent motion centered about a statistical mean value. This partitioning and ac- 
companying ensemble averaging leads to a set of Reynolds-averaged Navier-Stokes 
(RANS) equations. Although this process eliminates the need to completely resolve 
the turbulent motion, its drawback is that unknown single-point, high-order corre- 
lations appear in both the mean and turbulent transport equations. The need to 
model these high-order correlations is the well-known closure problem. The RANS 
method is a robust, easy to use, and cost effective means of computing both the 
mean flow as well as the turbulent stresses (velocity second-moments) and has been, 
overall, a useful flow prediction technology. 

The second approach is a filtering approach, and is most commonly known as 
the large eddy simulation (LES) method. Here the flow field variables are parti- 
tioned into resolved and unresolved scales. The original ideas were proposed over 
40 years ago by Smagorinsky (Smagorinsky 1963), and extended a few years later 
by Deardorff (Deardorff 1970). In these early works, the formal analogy between 
the unclosed subgrid scale stresses in LES and the unclosed Reynolds stresses in the 
RANS approach was exploited. The LES method is becoming a popular method for 
flow field predictions due to the rapid increase in computational power. While not 
as computationally demanding as DNS, proper implementation of LES - especially 
in the vicinity of solid boundaries - does require extensive computational resources. 
In addition, subtle issues associated with filter size and subgrid scale stress models 
are still being debated. What remains to be determined is whether LES can be 
used with sufficient confidence to provide the accuracy now being demanded of the 
more established RANS methods. It is clear, however, that as computational power 
increases the LES approach as well as variants to it (e.g. Geurts 2001) will become 
more prevalent. 

The increasing demand for more detailed information about the flow fields over 
complex aerodynamic configurations or in turbine engines, and the need to control 
such flows, has led to the need for methodologies capable of capturing some part 
of the instantaneous motion. These turbulent flows may still be statistically steady 
or stationary, but now only a portion of the turbulent motion needs to be modeled. 
For statistically unsteady flows, the need for alternatives to the RANS approach is 
more obvious since the foundation for RANS models is based on (spectral) equilib- 
rium hypotheses that may no longer be valid (e.g., Carpy & Manceau, 2006). The 
need for such predictions has led to the development of hybrid methods capable of 
having a RANS-type behavior in the vicinity of a solid boundary and an LES-type 
behavior away from the boundary. Such methods can also be interpreted within the 
usual LES context but with specially designed subgrid-scale (SGS) models capable 
of the dual behavior just described. The most often utilized approach to date is 
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the detached eddy simulation (DES) approach developed by Spalart and colleagues 
(Spalart et al. 1997, Spalart 2000). An insightful discussion of this topic, from an 
LES perspective, is given by Piomelli and Balaras (2002). Hybrid methods can be 
essentially classified in two categories: methods where the domain is decomposed 
into zones where LES is performed and zones where RANS is performed, with a 
sharp interface or, possibly, an small overlap region; or methods where the govern- 
ing set of equations is smoothly transitionning from a RANS behavior to a LES 
behavior, based on criteria updated during the computation. The models in the first 
category are often termed as zonal and the models in the second category as non 
zonal (Bunge, 2006), although this terminology is ambiguous since both are using 
different models in different zones. Both kinds of models will be a topic of interest 
and active research in the future. Currently, while such methods are often employed, 
little research has been conducted into formulating a consistent mathematical the- 
ory where such blending of methodologies is employed. Although the presentation 
will be somewhat general at the outset to show how the various approaches, at 
least formally, have the same form-invariant equations, the distinction being in the 
closure of these equations, the focus of this review will be on the traditional RANS 
approaches. Indeed, even in the hybrid RANS/LES frame, RANS modeling remains 
a cornerstone: the hybrid models highly rely on a RANS model in the near-wall re- 
gion, such that the influence of pressure gradients on the boundary layer and the 
prediction of flow separation is completely determined by the performance of the 
RANS model. Moreover, using a RANS-type model as subgrid-scale model can be 
favourable to account for complex phenomena such as heat transfer, effects of ro- 
tation or stratification, etc. In industrial applications, where a cost reduction is 
requested, the use of coarse meshes is not compatible with the hypotheses of a 
simple inertial behavior of the subgrid-scales, and the use of RANS-like transport 
equations and/or complex constitutive relations is and will be an active research 
field (e.g., Chaouat and Schiestel, 2005). 

From a physical standpoint in RANS model development, the task is to charac- 
terize the turbulence. One obvious characterization is to correctly describe the evo- 
lution of representative turbulent velocity and length scales, an idea that originated 
over 60 years ago (Kolmogorov 1942). This is the physical motivation behind the 
development of turbulent closure models. There is a hierarchy of turbulence models 
currently available for solving aerodynamic flow problems. These classes include the 
differential Reynolds stress models (DRSM), the algebraic Reynolds stress models 
(ARSM), the nonlinear eddy viscosity models (NLEVM), and the linear eddy vis- 
cosity models (LEVM). While other alternatives may exist, these form the basis for 
the majority of methodologies used, and amongst these the linear eddy viscosity 
models are the most common. 

The most commonly used linear eddy viscosity models in aeronautics are the 
Spalart and Allmaras (SA) model (Spalart & Allmaras 1994) and the shear stress 
transport (SST) model (Menter 1994), but in general purpose codes, different forms 
of high-Reynolds number K — e models, such as the standard K — e model (Laun- 
der & Spalding, 1974), the RNG K — e model (Yakhot & Orszag, 1986) and the 
K — e model of Shih et al. (1995a), often termed as “Realizable K — e” model, 
are still widely used. The main improvement of RNG and Shih et al. K — e model 
compared to the standard K — e model is the correction of the so-called stagnation 
point anomaly, through a modification of the e equation (RNG model) or of the 
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eddy-viscosity (Shih et al). The SA model is a one-equation model based primarily 
on empiricism and on dimensional analysis arguments. It is easily used with any 
type of grid: structured or unstructured, single block, or multiple blocks. The SA 
model has become popular among industrial users due to its ease of implementation 
and relatively inexpensive cost. The (SST) model has gained increasing favor due 
primarily to its robust formulation and improved performance for separated flows. 
It is a blend of the original K — uj formulation near walls and a I\ — e formulation 
in the outer region and in free shear flows. An important feature of the SST model 
is the modification to the definition of the eddy viscosity to account for the effect of 
the transport of the principal turbulent shear stress. More detail about the math- 
ematical form of these models can be found in the review by Gatski and Rumsey 
(2002). 

In aerodynamic flows where the turbulence is not confined to relatively narrow 
regions of the flow domain, such as on multi-element configurations at high angle of 
attack, wing-body junctions or turbine blades where separation or near-separation 
conditions exists, the demands on both accurate computational modeling of the 
physical problem as well as the quality of the turbulence model increase significantly. 
In proportion, the role of the CFD practitioner becomes more critical since poor 
modeling of the physical problem and/or poor choice of turbulence model can lead 
to correspondingly poor predictions. 

These increased demands on turbulence model performance have also high- 
lighted inherent deficiencies in the models themselves. For example, the inability 
to consistently predict separation accurately or to properly account for pressure 
gradient effects have become apparent in trying to replicate multi-element airfoil 
flow fields at higher angles-of-attack (Rumsey et al. 1998; Rumsey & Gatski 2001). 
Other examples are the difficulties to correctly account for rotation effects (Jakirlic 
et al . , 2002) and wall heat transfer (Thielen et al., 2005), which are both relevant 
to turbomachinery. A new deficiency has also been identified that is not directly 
related to turbulence model prediction itself, but to the application of turbulence 
models in predicting aerodynamic flow fields. This deficiency is related to the tur- 
bulence model being properly sensitized to transition onset location. While turbu- 
lence models are calibrated for fully turbulent flow predictions, aerodynamic flows 
include regions where the flow may be laminar and transitioning to fully turbulent. 
In attempting to predict and replicate such flows, a single system of equations is 
used throughout the computational domain and these equations are expected to 
predict the entire flow field which can include regions of laminar flow as well as 
fully turbulent flow. Some of these issues will be discussed further here. 

2. The Filtering Process and Governing Equations 

As detailed in the last section, it is currently necessary for numerical calculation 
of practical engineering turbulent flow fields to solve a set of equations for flow 
variables that represent the motion of a limited spectral range of scales. This de- 
scription holds true for LES, RANS formulations and any of the newly developed 
hybrid or composite methodologies currently being proposed. As such, the equa- 
tions describing the filtered motions in any of these formulations are form-invariant. 
They obviously differ with respect to the flow field motions being described, and 
this is predicated on how the higher-order correlations are parameterized. 
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As is customary in all these formations, the flow variable / is decomposed into 
a filtered part, /, and a sub-filtered part, /', given as 

/ = / + /'• ( 2 - 1 ) 

Generally, the filtering process can be defined as a subset of the general operation 
(e.g. Sagaut 2006) 

f(x,t) = g*f = J Q(x - x',t - t') f(x',t') dx' At' . (2.2) 


Different forms for the convolution kernel can be associated with the various solution 
methodologies. For the Reynolds-averaged formulation, for example, stationarity is 
usually assumed and a long time average then corresponds to an ensemble average. 
The filter function in Eq. (2.2) is given by 

£(x, t) = £ r (x, t) = G(x) G T (t) = <5(x) ^H(T - t) H{i) , (2.3) 

such that 

7 r (x, ^ /(x, t')dt' . (2.4) 

In such flows, the entire spectral range of scales is modeled so the sub-filtered part 
f is a fluctuating quantity whose average is zero /' = 0, and the mean quantity / 
can be extracted from 

1 f T 

E{f(x,t)}= Inn / t (x,T)= Inn — / f(x,t)dt, (2.5) 

T — kxj T — >oo 1 J q 


Note that, since limr^oo Gr(t) = 0, Reynolds averaging cannot be expressed as a 
convolution filter, although the Reynolds averaged function is the limit of a series 
of functions obtained by convolution, as shown by Eq. (2.5). As a practical matter, 
Reynolds averaging is evaluated as long-time averaging, where “long-time” is not 
infinity but rather a time sufficiently large compared to the turbulence time scale. 
Therefore, Reynolds averaging can be considered, within any predefined accuracy, 
as the convolution filter Gt with a sufficiently large T. A useful property of this 
filtering is that the average of the product of two quantities is f g = f g + f'g'. 
This allows for the easy extraction of stationary correlation data from numerical 
simulations where the instantaneous values of / and g are computed. 

For a flow statistically periodic in time (cyclo-stationary), it is customary to use 
phase averaging, corresponding to a filter function given by 


£r(x,f) 


G(x) G r (i) 


<*(x) 


1 


N 


lim 

N — »oo N + 1 


Y^Sit + nT) 


n — 0 


( 2 . 6 ) 

(2.7) 


where T is the period of the cycle. In that case, the phase average of the sub-filtered 
part /' is zero, < f >= f = 0, and the filtered, or phase-averaged quantity < / > 
can be extracted from 


< /(x,f) >= /(x,f) 


Gt * f 


lim 

N — xx) 


l 

N + 1 


N 


Y.n^+ n n 

n—0 


(2.8) 
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For the most part, the large eddy simulation methodology has been based on 
spatial filtering. The convolution filters most commonly employed are the top-hat 
filter, the Gaussian filter and the spectral cutoff filter. In numerical algorithms 
requiring structured meshes applicable to multi-dimensional flows, such spatial fil- 
tering can become problematic since anisotropies in the filtering process need to be 
introduced. For unstructured grid schemes, the anisotropy issue is compounded. An 
alternative approach that helps to alleviate this problem is the Eulerian temporal 
filtering approach. (It is recognized that the problem cannot be completely solved 
since an inherent filtering occurs even on the spatial scales with a temporal filtering 
process.) This procedure has received little attention in the past although some re- 
cent work (Pruett 2000; Pruett et al. 2003) suggests it can be an effective alternative 
to the spatial filtering process. Causal time domain filters can be constructed that 
are analogues to the spatial filters. A simple example of causal filter is the top-hat 
filter given by Eq. (2.3) where T is the temporal filter width. Analogous to the spa- 
tial filtering approach, only a portion of the spectral range of scales is modeled with 
sub-filtered part /'. The filtered quantity / can then be extracted from Eq. (2.4). It 
can thus be seen that within the realm of temporal filtering, it is possible to develop 
a more rigorous linkage between the large eddy and Reynolds-averaged approaches 
(Pruett et al. 2003), since the Reynolds-averaged function E{f} is the limit of the 
temporally-filtered function fj- when the temporal filter width T goes to infinity. In 
the frame of spatial filtering, such a linkage can only be established in homogeneous 
flows. 

The aerodynamic flows of current relevance are compressible and the relevant 
governing equations, for whichever solution methodology is used, should be cast in 
the appropriate form. This is becoming even more important as the interest in high- 
speed flow fields increases. It should be noted, however, that having a turbulent flow 
compressible does not mean the turbulence dynamics is necessarily compressible. 
For attached flows, mean or filtered flow fields over most of the supersonic regime 
(M < 5) are described reasonably well by using incompressible dynamics under 
adiabatic conditions and in the absence of shocks. 

The filtering processes described above yields governing equations that are all 
formally equivalent. They differ, not in form, but in the models or parameterizations 
needed to close the equations. The formal similarity of these equations, irrespective 
of the filtering process invoked to derive them, is the foundation upon which hybrid 
or composite methods can be constructed. However, the modeling or parameteri- 
zations of the various terms requiring closure is a difficult and challenging problem 
since the scales of motion being represented by each filtering procedure can be dif- 
ferent. In addition, for compressible flows it has been found that a rewriting of the 
equations using mass- weighted (Reynolds 1895), or Favre variables (Favre 1965) is 
advantageous, since the equations take a more compact form and are structurally 
similar to their incompressible counterpart. 

For a dependent variable /, the Favre mean is defined as 


/ = 



(2.9) 


The instantaneous value / can then be decomposed into either the usual Reynolds 
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averaged variables or the Favre-averaged variables 

/ = / + /' = / + /"• ( 2 - 10 ) 

As might be expected, an extensive list of relations exist between the Reynolds- 
averaged variables and the Favre variables (Barre et al. 2001). 

The equations for the mean density p and mean momentum pfq are given by 

% + s- <7ii > , = 0 - (211) 


_Diii 

P ~Dt 


d{pui) 

dt 


d __ dp dan 


d{pTij) 

dxi 


and the mean viscous stress tensor <j, ? is 


&ij — 2/i I Sij 0 Skk^ij — 2/i ( Sij 0 Skk 


(2.12) 


(2.13) 


where JJ is the mean molecular viscosity, and Ty = it" it " is the Favre-averaged 
velocity correlation tensor. Equation (2.13) neglects contributions from pi , and 
assumes that iq « iq. This assumed equality between the average velocities implies 
that the average fluctuating velocity u" is small since iq — u i = it". 

The equation for the mean total energy pE is 


where 


and 


^tst + £- fo 6 ) = fa + pE "<) ■ 


E = c v T - 


UiUi , u"u" 


H = E+ 

P 

, dT - dT 
q, = -fcr ^ — - - k T a \ — 
J dxj dxj 


(2.14) 

(2.15) 

(2.16) 

(2.17) 


pE"u'\ = CppUjT" + Ui (pTij - (Tij) + 


pu'iU'iU " 


- CTi 


U i 


(2.18) 


In the approximation of (2.17), fluctuations in the thermal conductivity are ne- 
glected, and the Favre-averaged and Reynolds-averaged mean temperatures are 
taken as approximately equal. The equation of state in mean variables is 


p = pRT, 

or in terms of the mean total energy ~pE 


(2.19) 


P={ 7-1) 


pE — -p (it 2 + v 2 + w 2 ) — pK 


(2.20) 
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where 7 is the ratio of specific heats ( c p /c v ), and R is the gas constant. The presence 
of the turbulent kinetic energy term K . 

pK = p^=p^-, (2.21) 

suggests a strong coupling between the mean equations and the turbulent transport 
equations. 


3. Closure Strategies for RANS 

The closure problem arises at all levels of the statistical moment equations. At 
the mean flow level, the RANS momentum equations Eq. (2.12) require closure 
through a specification of T,y, and the mean energy equation Eq. (2.14) requires the 
closure of the turbulent heat flux pc p it"T". At the second- moment level, closure is 
required for the velocity-pressure gradient correlation, tensor dissipation rate, the 
velocity-triple moments, and the turbulent mass flux it" = —p'u'/~p. Closure for 
these terms is achieved either through the solution of partial differential transport 
equations or by assuming a tensor polynomial expansion for the unknown corre- 
lation in terms of basis tensors formed from the independent tensors on which it 
depends (Gatski 2004). It is not possible in the space allotted here to discuss the 
modeling of all these unknown higher-order correlations. Only the correlations re- 
lated to the velocity-pressure gradient correlation and the tensor dissipation rate 
will be discussed. These correlations are relevant to the solution of both incom- 
pressible and compressible flows and, in the case of the tensor dissipation rate, it 
has been found that the solenoidal (incompressible) dissipation rate need only be 
accounted for (e.g., Kreuzinger et al. 2006). It is safe to assume that the newer 
hybrid methodologies have not reached a level of maturity where detailed modeling 
issues associated with many of these correlations have been addressed. 

(a) Levels of Turbulence Closure 

It is often desirable to work with the turbulent anisotropy tensors or scaled 
scalar flux vectors in analyzing turbulent flows. For the Reynolds stress tensor, the 
anisotropy tensor bij = ( Tij—2K8ij/3)/(2K ), where K = tu/2 , can be defined, and 

for the heat flux vector, for example, the scaled scalar flux vector u”T" /(KT' /2 ) 
can be defined. 

The differential Reynolds stress model requires the solution of six partial differ- 
ential equations for the Reynolds stress components and a transport equation for a 
variable from which a length scale can be obtained. For the most part, this scale re- 
lated equation is the turbulent kinetic energy dissipation rate equation. In addition, 
various scalar flux equations, each having three components, are also required. As- 
sociated in general with the scalar flux equation is the need for a scalar variance and 
corresponding scale equation for the scalar variance. In general, the full differential 
form at the second-moment level can be computationally challenging. 

Subsets of this full differential set are often desired and can be obtained. Implicit 
algebraic Reynolds stress and scalar flux equations can be obtained directly from 
the differential forms by invoking what are termed weak equilibrium assumptions. 
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Applicable models are developed through explicit polynomial expansions (represen- 
tations) of the anisotropy tensor and scaled scalar flux vector. The proper choice of 
basis is dependent on the functional dependencies associated with the anisotropy 
tensor and scalar flux vector. Note that the equations for the Reynolds stress 
anisotropy bjj are only frame-indifferent under the Galilean group of transforma- 
tions. Under the more general Euclidean group, the equations exhibit a dependence 
on system rotations and translational accelerations. Transport equations are still 
required since the algebraic models developed simply yield an expression for the 
components of the Reynolds stress anisotropy tensor and the scalar (heat) flux 
vector. These transport equations are associated with the turbulent kinetic energy 
and the energy dissipation rate (or dissipation rate per kinetic energy), and scalar 
variance and scalar variance length scale. 

What are termed nonlinear eddy viscosity models form a link between what one 
would call the higher-order models - the differential and algebraic Reynolds stress 
models, and the lower-order linear eddy viscosity models. The nonlinear eddy vis- 
cosity models describe the turbulent Reynolds stress field by a polynomial expansion 
similar to the algebraic stress model; however, the expansion coefficients are deter- 
mined from calibrations with experimental or numerical data, and on some physical 
consistency constraints (e.g., Shih et al. 1995b; Craft et al. 1996). This contrasts 
with the algebraic stress models where the expansion coefficients are extracted di- 
rectly from the Reynolds stress transport equations. An analogous situation holds 
for the scaled scalar flux vector. Once again, transport equations are required as in 
the case of the algebraic anisotropy tensor and algebraic scalar flux vector repre- 
sentations. 

At the linear eddy viscosity and diffusivity level of closure, the mean momentum 
and temperature equations are closed by using a Boussinesq-type approximation 
between the turbulent Reynolds stress and the mean strain rate tensors, and scaled 
scalar flux and gradient of the mean scalar property. Although computational re- 
sources have become less of a prohibitive factor in choosing turbulence models, 
many practitioners have continued to rely on the lower-order eddy viscosity and 
diffusivity models. Much of this is due to robustness issues associated with the 
higher-order models as well as to lack of familiarity with the formulations. 


( b ) Improved High-Order Correlation Closures 

The main focus of advanced RANS turbulence modeling research for aerody- 
namic flows is on correlations associated with the turbulent velocity field in isother- 
mal flows. This does not suggest that modeling issues associated with the scalar 
fluxes (e.g. heat and mass) are unimportant, rather it is a reflection on the limited 
scope of advanced modeling activity currently being supported. For the turbulent 
velocity field, the differential Reynolds stress transport equations are the starting 
point for model development. The two unknown correlations that currently receive 
the most attention are the velocity-pressure-gradient and tensor dissipation rate 
correlations. These terms are important because they represent the mechanisms 
responsible for the redistribution of normal and shear stress throughout the entire 
boundary layer flow including the near-wall region. 
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(i) Pressure- Strain Rate Correlation 

In the case of the velocity-pressure-gradient correlation, it is customarily rewrit- 
ten in terms of the pressure-strain rate correlation and the pressure-velocity corre- 
lation as 


+ U '0£ = ( M + ^1) + JL (pV ) + — (v’u 1 ) ( 3 1) 

; dxj + j dxi P \dxj + dxi J + dxj [P + dxi ^ ^ ( ’ 


— Fli- 


Df, 


Initially, this partitioning was intended to isolate the effects of the pressure trans- 
port terms which predominantly contribute to the dynamics near solid boundaries. 
The pressure-strain rate correlation term Hy is then decomposed into a rapid and 
slow part. The rapid part is modeled by assuming dynamic equilibrium conditions 
in both physical and spectral space, and the slow part is modeled by assuming either 
a linear or nonlinear relationship with the Reynolds stress anisotropy tensor. The 
wall proximity corrections required will be discussed separately below. As might 
be expected, an extensive literature base deals with the full range of issues asso- 
ciated with developing specific models (Hanjalic & Jakirlic 2002). Almost all the 
research efforts have focused on the pressure-strain rate correlation since this term 
is of the same order as the production term and acts as a redistribution between 
the Reynolds stress components. In compressible flows without shocks and close 
to adiabatic conditions, the partitioning in Eq. (3.1) is applicable and for the most 
part the pressure-velocity term is neglected. In compressible flows with shocks, the 
pressure- velocity correlation has a non-negligible effect on the dynamics even away 
from solid-boundaries and should be accounted for. 

Away from solid boundaries, the pressure-strain rate correlation also acts to 
diminish the difference between the normal stress components. In the vicinity of 
the solid-boundaries, its action, along with pressure-velocity correlation, increase 
the anisotropy of the stress field. This action along with that of the tensor dissipa- 
tion rate, is to enforce the two-component limit on the Reynolds stress tensor. In 
aerodynamic flows of interest, these inhomogeneous effects can be important, so it 
is clear why such terms have received so much attention. 

In all of the commonly used second-moment closure models, the pressure-strain 
rate correlation 11^ is modeled away from solid-boundaries in the general form 
(Chou 1945; Lumley 1978) as 


n v = eAij (b) + KMijki (b) , (3.2) 

where A t] (b) and A4iju{b) are tensor functions of the anisotropy tensor and are 
related to integrals over the flow volume derived from a pressure Poisson equation 
for incompressible flows or a convective pressure wave equation for compressible 
flows, and e is the isotropic dissipation rate. The isotropic dissipation rate is by 
the second law of thermodynamics always positive, and this Reynolds-averaged 
turbulent dissipation rate can be written as 

F?) (3.3) 


pe — 2/i ( [s':s'] — 
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where s' = sE = (uM +u' i )/2, { } is the trace, and [ : ] is the trace of the product. 
The second term on the right is associated with fluctuating dilatation and can be 
neglected (Pantano & Sarkar 2002; Pirozzoli et al. 2004), and the first term can be 
related to the fluctuating enstrophy (e.g. Sarkar et cd.1991). 

In Eq. (3.2), the Aij(b) term is usually associated with the “slow” relaxation of 
the turbulence toward isotropy, and the Mijki{b) term is usually associated with 
the “rapid” response of the turbulence to imposed mean velocity gradients. This 
partitioning has its origins in incompressible flows where the turbulent pressure 
held p' is split into slow p'^ and rapid pd R ) parts. The pE s l part is the solution 
of a Poisson equation that only involves gradients of the turbulent velocity held; 
the pd R ) part is the solution of a Poisson equation involving the mean velocity 
gradients. The terminology originated from the observation that the slow part will 
only adjust as the turbulence itself adjusts, but the rapid part will adjust instantly 
through the mean velocity gradient. 

For statistically stationary turbulence, the starting point for the model devel- 
opment of the rapid part can be treated with the transform pair 

IIij{ r) = p'(x + r, t ) t) + -j^u" (*, t)^ , (3.4a) 

and its equivalent Fourier transform, 

JTjj(r) = y J d 3 xj d 3 k'd 3 kV k ' (x+r V k " x 

x (p(k') (k'f Ui( k") + k'luA k"))> , (3.4b) 

where r is the two-point separation distance (r — > 0 for single point closures), u 
and p are the transformed huctuating velocity and pressure helds. From this equa- 
tion, it is now possible to obtain an expression for the rapid part of the pressure- 
strain rate correlation once an appropriate representation for the energy spectrum 
tensor is found. (The slow part of the pressure-strain rate correlation, represent- 
ing turbulence-turbulence interactions, yields triple- velocity terms. Models for this 
“slow” contribution are derived in an alternative manner by simply relating it to 
the Reynolds stress anisotropy.) 

For incompressible flows, the fourth order tensor M.ijki (b) given in Eq. (3.2) is 
related to the energy spectrum tensor through the simple relation 

M ijpq (b) = J ^Eiji b,k)d 3 k , (3.5) 


where the energy spectrum tensor E l:] (b, k) can be represented by a polynomial ex- 
pansion involving both the Reynolds stress anisotropy and the wavenumber vector. 
The simplest representation that can be used is given by four-term representation 


Eij (k, b) 


m 

Airk 2 


( kikj \ E a (k) fk n k m \ ( k l k j \ 

(*« - — J + (— H ( j « + “P" ) 

, E a (k) r /, k n kj kik n \ 


(3.6) 


where E{k) is the isotropic energy spectral density and E a (k) is the anisotropic 
energy spectral density. 
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With the functional dependency outlined in Eq. (3.2), the proper representation 
for II would be generated from the integrity basis given by the invariant combina- 
tions of b, S, and W (e.g. Gatski 2004). Obviously, any representation composed 
of the full integrity basis would be unmanageably large. Even though the energy 
spectrum tensor given in Eq. (3.6) is only a linear function of the anisotropy tensor 
b, many representations for the rapid part of the pressure-strain rate correlation 
have used bases with quadratic and cubic terms in b (cf. Sjogren & Johansson 
2000) although they only retain a linear dependence on the mean velocity gradient 
field. This linear dependence on the mean velocity gradient is consistent with the 
functional form for the pressure-strain rate correlation given in Eq. (3.2); however, 
appearance of the higher-order powers of the anisotropy tensor is inconsistent with 
the linear dependence on the anisotropy tensor assumed for the energy spectrum 
tensor. 


(ii) Dissipation Rate Tensor 


The other unknown correlation that appears in the differential Reynolds stress 
equations and that continues to be of modeling interest is the tensor dissipation 
rate ey. Even though it is possible to derive a transport equation for J>e v j from the 
fluctuating momentum equation by taking the moment (cf. Speziale 1991) 


2/i 


MJL (kt >\ du 'i d 

dx k dx k Uj + dx k dx k 


(A/Y.) 


= 0 , 


(3.7) 


where Af is the Navier-Stokes differential operator, the final form of this equation is 
quite complex and contains terms where very little, if any, information is available 
for closure modeling. 

It has been assumed, that away from solid boundaries, the tensor dissipation 
rate becomes isotropic (e,j = 2/3 e <5,j); although, there have been suggestions that 
such assumption may not be generally true (Durbin & Speziale 1991). Thus, early 
attempts at modeling the deviatoric part of the dissipation rate tensor focused on a 
coupling with the Reynolds stress anisotropy (Hanjalic & Launder 1976; Hallback 
et al. 1990). 

The quantity ~pe represents more than a transformation of kinetic energy to 
internal energy since it can be further decomposed into 


pe = n 


. dut du'j 
dx k dx k 


d 
dx k 


du' 


dxi 


pe’ 


(3.8) 


The second term, ~pe' , is the divergence of a flux, and as such is purely a transport 
term since it redistributes energy from one region to the other, but does not con- 
tribute to the global time evolution of kinetic energy. This term has been evaluated 
by Bradshaw & Perot (1993) in a channel flow and represents less than 2% of the 
total. 

The energy dissipation rate also plays a pivotal role in the dynamics of the 
turbulent kinetic energy. It can be extracted from the viscous term, 


-pT^jiur 


<9 2 u' 

dxjdxj 


(3.9) 
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which is not necessarily positive and requires modeling, that appears in the trans- 
port equation for the turbulent kinetic energy. It is customary to decompose this 
term into 

iiT '' = p a|£"“' i£ - <3 - 10) 

s. v ✓ 

~pD v 

where the first term, ~pD^ is called turbulent diffusion. This decomposition has two 
very desirable properties: pD^ does not require modeling, and pe is always positive 
which facilitates its modeling. 

The variable e is often called the homogeneous dissipation, since pT' 3 = ~ps = ~pe 
in homogeneous flows, and e' = e — e is sometimes called the inhomogeneous dissipa- 
tion. It is worth emphasizing that the decomposition of homogeneous/inhomogeneous 
dissipation is completely artificial, since there is an infinity of possible decompo- 
sitions: for instance, any linear combination of the form e + can be called 
homogeneous dissipation, since it goes to e in homogeneous flows. In the differential 
Reynolds stress formulation, where the isotropic assumption for the dissipation rate 
tensor is used, this isotropic dissipation rate is obtained from a modeled transport 
equation (Hanjalic & Launder 1976). 

There have been recent attempts to account for dissipation rate anisotropy 
through a tensor dissipation rate equation (e.g., Oberlack 1997; Speziale & Gatski 
1997, Jakirlic & Hanjalic 2002). While not explicitly focused on developing closures 
to account for wall proximity effects, these studies used tensor representations to 
extract either an explicit algebraic dissipation rate model or a differential model 
for the tensor dissipation. Unfortunately, extensive validation studies are difficult 
due to the lack of data available for a quantity such as . Numerical simulations 
provide the best source of data for such quantities but are unfortunately limited to 
simpler flow geometries. 

It is now recognized that effects of anisotropies in the turbulence statistics should 
be accounted for in some way for many flows of practical interest. Obviously, the 
full differential models do this, but at a higher computational overhead than the 
linear eddy viscosity models. One of the main reasons that algebraic vector and 
tensor polynomial representations have become popular is that such anisotropies 
can be accounted for at only a small increase in computational cost over the linear 
eddy viscosity models. The details of this representation procedure have been put 
forth previously (e.g. Gatski 2004, Gatski & Wallin 2004) and it is not necessary 
to repeat them here. 


(c) Modelling of Near- Wall Turbulence 

While the discussion up to this point has been focused on mathematical repre- 
sentations of the rapid pressure-strain rate correlation in the case of homogeneous 
flows, the practical application of such models to inhomogeneous flows requires 
accounting for the effects of solid boundaries. This inherently brings into consid- 
eration the proper accounting of turbulence anisotropy into the modeling process. 
Although this effort has sometimes been called low Reynolds number modeling, 
it is more properly called near-wall modeling to distinguish it from the separate 
problem of transition prediction discussed in the next section. 
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Both the pressure-strain rate correlation and the anisotropic dissipation rate 
that appear in the transport equation require modification, and the pressure- 
velocity correlation cannot be neglected anymore. Along with these parameteriza- 
tions, which focus on the stress transport equations, is the behavior of the dissipa- 
tion rate equation itself in the vicinity of the wall. Overall, the issue of near-wall 
modeling is at least as complex as the issue of developing the high-Reynolds-number 
models, and unfortunately less precise. 

Several attempts have been made to develop near-wall closure corrections for 
the Reynolds stress transport equations and the dissipation rate equation. For the 
turbulent stress transport equations, attempts have almost exclusively focused on 
extending the high Reynolds number pressure-strain models. This approach has 
been taken by Launder and co-workers who have developed a methodology that 
enforces the two-component limit of the turbulent Reynolds-stress field as the solid 
boundary is approached (Craft & Launder 2001). For the dissipation rate equation, 
additional terms and modified coefficients have been proposed to account for the 
anisotropic near-wall effects and the correct limiting behavior at the wall (e.g. 
Speziale & Gatski 1997, Jakirlic & Hanjalic 2002). 

A major focus over the last decade has been on what is termed the “elliptic 
relaxation method” (Durbin 1991, 1993; Manceau & Hanjalic 2000; Manceau et al. 
2001; Manceau et al. 2002; Laurence et al. 2005). Based on a theoretical analysis of 
the influence of the blocking of the wall, this approach introduces a tensor function 
representing the combined effects of a near-wall velocity-pressure gradient corre- 
lation and anisotropic dissipation rate. The tensor function asymptotes to a high 
Reynolds number form away from solid boundaries through an elliptic relaxation 
equation. 

It is worthwhile to briefly outline the formulation. The transport equation for 
Tij can be written in the form 


Drij_ 

dt. 


- dT ij D 
‘ Uk 0 — Pij 

dx k 


- £i 


D l 

v 


d 2 i 


dx k dx k 


(3.11) 


where (pij is not traceless (cf. 77y ) , since it includes pressure diffusion (see Eq. 3.1). 
In the elliptic relaxation method, the variation of the dissipation rate anisotropy dij 
( dij = kij — 2e5ij/3)/(2e)) as the wall is approached is accomplished by a relaxation 
to its wall value, which is assumed to be equal to bij . With this assumption, (pij — &ij 
in Eq. (3.11) can be written as 


(pij - £ij = eKfij - T -^e , (3.12) 

with the relaxation function /y defined by 

eKfij = (p^ - 2 e ( dij - b^) . (3.13) 


The original scaling of the relaxation function f tl was solely through the turbulent 
kinetic energy; however, Manceau et al. (2002) have shown that adding a dissipation 
rate factor e to the scaling (eKfij) eliminates an unwanted amplification effect 
inherent in the original scaling. The system is closed through the solution of a 
relaxation equation for fjj 

(1 - L 2 V 2 ) fij = P (77y + 2 ebij) , (3.14) 
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where H^ can in principle be any of the high Reynolds forms found in the literature. 
Although some authors have tested nonlinear models (Wizman et al. 1996), terms 
linear in bij are used in general. Moreover, despite the fact that pressure diffusion 
can readily be accounted for in Eq. (3.14), this term is usually considered negligible 
far from the wall, i.e. , in the right hand side. Thus, the elliptic relaxation equation 
is driven by the high Reynolds number form of the pressure-strain rate correlation 
FT and a contribution from the Reynolds stress anisotropy 2eb (away from the wall 
the dissipation rate is assumed to be isotropic dij = 0). 

It is interesting to contrast the predictive capabilities of previous near-wall cor- 
rections with the performance of the more rigorous elliptic relaxation approach. 
One of the first Reynolds stress models proposed to handle near-wall effects was 
developed by Hanjalic and Launder (1976). Figure 1 shows a comparison of results 
from the elliptic relaxation method (Manceau et al. 2002) and the Reynolds stress 
model (Hanjalic & Launder 1976) with direct simulation data at Re T = 180 for 
the Reynolds shear stress (in wall units). The results of the Speziale, Sarkar, and 
Gatski (SSG) model (Speziale et a!1991), without near-wall correction, is also in- 
cluded for comparison. At this low Reynolds number, the RSM approach based 



Figure 1. Comparison of wall proximity cor- 
rections to RANS closure for fully devel- 
oped turbulent channel flow. 



Figure 2. Isocontours of a 2 in a backstep 
flow. The solid lines correspond to the levels 
0.25, 0.5, 0.75 and 0.9. 


on a modification of the tensor dissipation rate and the introduction of an extra 
production term was not able to accurately predict the near-wall behavior of the 
Reynolds shear stress. The elliptic relaxation method, on the other hand, did an 
excellent job at this Reynolds number, and as shown in Manceau et al.( 2002) also 
made excellent predictions over a much larger Reynolds number range. 

However, Reynolds-stress models based on elliptic relaxation have not spread 
into industrial codes and applications, because of the huge additional numerical 
effort required compared to simpler closures (the vr-f model for instance, which 
includes elliptic relaxation in a linear eddy- viscosity framework, has become popular 
and is now available in several commercial codes). The reason does not only lie in 
the requirement of solving six additional differential equations Eq. (3.14), but also 
in the numerically problematic wall boundary conditions for the components / 22 , 


Article submitted to Royal Society 


16 


T. B. Gatski, C. L. Rumsey & R. Manceau 


fi 2 and / 23 : 


fa I =-20 — 
J L J \vj ,~2 


lim 

® 2 =o a;5 


(3.15) 


where 0:2 is the wall-normal direction. (The boundary conditions for the other com- 
ponents are: f\\ = / 3 3 = — 5/22 and /i 3 = 0.) One of the difficulties is that these 
boundary conditions must be applied in the local frame determined by the orien- 
tation of the wall, which couples the Reynolds stress equations via the boundary 
conditions when the wall is not aligned with the global reference frame. 

Manceau and Hanjalic (2002) have proposed the elliptic blending Reynolds- 
stress model (EB-RSM), which tries to preserve all the desirable properties of the 
elliptic relaxation model, in particular, those properties related to the representation 
of the blocking effect of the wall; while decreasing the number of equations and the 
numerical stiffness of the boundary conditions. In the most recent version of the 
model (Manceau 2005), (pij — Sij is modeled as 


<t>ij - £ ij = (1 - « 2 ) (4>7j - £ ij) + « 2 - £ ij ) . (3.16) 


The purpose of the introduction of the near- wall form <(>)" — is the reproduction 
of the asymptotic behavior of (f>ij — e^, a role that is played by the boundary 
conditions Eq. (3.15) in the elliptic relaxation model. Here, the near- wall form of 
the dissipation tensor is chosen as efj = T-ijs/k, but any other choice is possible: 
the important point being that to satisfy the behavior of <f>fj — e^-, the model for 
<f>fj must be adapted to the model for In that case, <j>ij must take the form 


/ — u ,2 /2 uV 0 \ 

4>ij ~ 4"Tj ~ n K £ fij I-UJ = — ^ 7 7 I u'v' v' 2 v'w' I (3-17) 

V 0 v’w’ - v ' 2 /2 J 

Contrary to boundary conditions that can be imposed in a local frame linked to 
the wall where they are applied, Eq. (3.16) is applied inside the flow domain, and 
requires information about the wall distance and the orientation of the wall. The 
wall distance is felt implicitly by the blending function a , which goes from 0 at the 
wall to 1 far from the wall, since it is obtained via an elliptic relaxation equation, 
very similar to the ones used in the elliptic relaxation model 

a — L 2 V 2 a = 1 , (3.18) 


with the boundary condition a = 0 at the wall. In a semi-infinite domain bounded 
by a flat plate, the analytical solution of this equation would be 

a = 1 — exp (^j^j (3.19) 

if L was taken as a constant. Therefore, a is a function of the ratio y / L=Wall 
Distance / Integral Length Scale that the model feels. Figure 2 shows isocontours of 
a 2 around the step corner in the case of a backstep flow at Re = 37500. It can 
be seen that the width of the region over which the near-wall model is active is 
strongly dependent on the local flow conditions, via the length scale L, which is 
much larger in the recirculation region than in the incoming boundary layer. 
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The second important issue is the orientation of the wall. It is crucial that the 
model is valid regardless of the relative orientations of the wall and the reference 
frame, i.e. , the near- wall model must be objective. Indeed, since the (f>ij tensor is 
objective, it is crucial to ensure that <t>™j is as well. This requirement is automatically 
ensured if (/>)" is made a function of objective quantities. In order to correctly enforce 
the asymptotic behavior Eq. (3.17), it is seen that the model cannot be proportional 
to Tij , but must rather be able to distinguish between the components, depending 
on the orientation of the wall. Therefore, an objective tensor must be defined that 
provides information about the orientation of the wall. First, it can be noted that 
the vector n = Va is normal to the wall in its vicinity, since the wall is the 
isocontour a = 0, and is objective, since it is the gradient of an objective scalar. 
Thus, the tensor 

N = n ® n = (■ UiUj ) (3.20) 

is also objective, and a simple way to reproduce Eq. (3.17) is to use the model 


b w — — 

b K 


tN + Nt- — [t:N] (N + I) 


(3.21) 


This near-wall form is then transitioned to the high Reynolds number form 

fiij ~ £ ij = n ij — ■ (3.22) 

Here, the SSG model (Speziale et a/.1991) is used for II, although any high Reynolds 
number model could be used. 

Thus, with such a blending approach, the distance to the wall and its orientation 
are felt by the model through a and the tensor N. Since no geometrical information 
is explicitly introduced in the equations, the model is easy to use in complex ge- 
ometries. Moreover, the number of equations necessary to enforce correct near-wall 
behaviors is reduced from 6 to 1, as compared to the elliptic relaxation model, and 
the wall boundary conditions are much simplified, which is crucial from a numerical 
stability point of view. It is worth pointing out that the model is able to reproduce 
the near-wall anisotropy of turbulence as well as, or better than the full elliptic 
relaxation model, as shown in Fig. 3. 

Another focus of near- wall research, that has and continues to be of interest, is 
the modeling of the near- wall anisotropy dij of the dissipation tensor e t j . Contrary 
to the elliptic relaxation/blending strategies, discussed above, in which the differ- 
ence (f>ij — £ij is modeled as a whole, many authors have investigated separately the 
dissipation tensor. Provided that all the models are based on an isotropic dissipa- 
tion far away from solid boundaries, the modeling challenge consists of providing a 
smooth transition towards an anisotropic near-wall formulation. 

The decomposition Eq. (3.10) between viscous diffusion and homogeneous dissi- 
pation discussed previously is not without modeling implications for the dissipation 
tensor. Indeed, all the models are based on an algebraic formulation of the dissi- 
pation tensor, which ensures a smooth transition between a high Reynolds form, 
applied in regions far away from solid boundaries, and a near-wall form, derived in 
order to reproduce the asymptotic behavior of the dissipation tensor in the vicin- 
ity of the wall. In the standard decomposition, where the homogeneous dissipation 
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Figure 3. Comparison of Reynolds stress 
anisotropies predicted by the Elliptic Re- 
laxation Model (13 equations) and the El- 
liptic Blending Model (8 equation) 


Figure 4. Comparison of the model f e bij for 
the dissipation anisotropy with two possi- 
ble definition dij and d*j. A priori test in 
a channel flow at Re r = 590 (Moser et al. 
1999). 


tensor e-ij is defined by 

„ du'jdu'j 

^ dxjdxj ^ dxk dxk 

pD <*. P e » 

r tj 


(3.23) 


the asymptotic behavior of the anisotropy of the dissipation dij is difficult to relate 
to the anisotropy of the Reynolds stress 6 ,, . Indeed, for a wall located in X 2 = 0, 
the components of d l: j and b l; j are asymptotically related by 


d^ 


bn 2 6 12 &13 \ 

2 612 4 622 2 623 J 

bi 3 2 623 633 / 


(3.24) 


which shows that the two tensors are not proportional. Therefore, the standard way 
of modeling the dissipation anisotropy 


dij — f e 6 


13 > 


(3.25) 


where f e is any function approaching unity at the wall and going to zero far away 
from the wall, does not provide the correct limiting anisotropy for the components 
12, 23 and 22. 

The analysis of the transport equations of the two-point correlation tensor 
u i A u 'j B (Jovanovic et al. 1995) suggests that the decomposition of viscous diffu- 
sion/homogeneous dissipation should be written as 


p = 1 _ d 2 u'ju'j _ du' i du'j 1 _ dPuty 

^ ^ 2 ^ dxjdxj ^ dxk dxk 2 ^ dxjdxj 

s. ' ,✓ s. 


(3.26) 


Jakirlic & Hanjalic (2002) have pointed out that with this decomposition, it is much 
more justified to assume proportionality between the anisotropies 


djj — febij , 


(3.27) 
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since the asymptotic analysis leads to 

bn b \2 &13 \ 

612 2 &22 &23 I (3.28) 

6 13 &23 &33 / 

although the 22 component is still not exact. Figure 4 shows a comparison of the 
model f e bij with the dissipation anisotropy obtained from the DNS of Moser et al. 
(1999), based either on the standard decomposition Eq. (3.23) or on the decompo- 
sition Eq. (3.26) proposed by Jakirlic & Hanjalic. The blending function is modeled 
as f s = 1 — \/~AE 2 , where A and E are the flatness parameters associated with the 
Reynolds stress and the dissipation, respectively (Jakirlic & Hanjalic 2002). It can 
be seen that the model better represents d \ jj than dij , although, quite surprisingly, 
the improvement is especially visible on the components 11 and 33. 

These types of approaches, focusing either on terms involving the pressure fluc- 
tuations or on the tensor dissipation behavior, have yielded encouraging results in 
an effort to more accurately account for near-wall behavior, and improve the ro- 
bustness of such numerical predictions. As with the other closures discussed in this 
paper, only extensive validation tests will confirm their ability to handle a wide 
range of flow fields. 



4. Some Current Challenges 

The challenge problems to be briefly discussed here are by no means complete. Ar- 
guably, other application areas are of current high interest, including buffet onset 
and compressible mixing. However, the two rather broad areas discussed here have, 
in our view, been pacing items in the migration from traditional RANS methodolo- 
gies to hybrid RANS-LES methods and ultimately to large-scale simulation method- 
ologies. 


(a) Accounting for Transition and Re-laminarization 

Many aerodynamic flow fields of interest, where single-point closures based on 
the Reynolds-averaged Navier-Stokes methodology are applicable, describe regions 
where the flow is transitioning from laminar to turbulent or re-laminarizing from 
the turbulent state. RANS models have routinely been applied to such regions and 
the accuracy of the flow predictions in these “non-turbulent” regions have varied 
significantly depending on the models. Obviously, the RANS turbulence models 
were never calibrated to predict such transitioning regions; however, the current 
demands on RANS predictions of aerodynamic flows requires that such models 
predict with some degree of accuracy, the flow in these transitioning regions. 

At the outset, it is necessary to recognize that what is required of a RANS model 
is not the ability to predict the detailed transition dynamics per se, but rather that 
the RANS model be properly sensitized to predict the statistical characteristics of 
the flow in these transitioning regions. The basis for many of the models being used 
today originated with the ideas of Emmons (1951) (see also Dhawan & Narasimha, 
1958). The desire to account for transitional effects in RANS models has been given 
renewed emphasis and has been ongoing for some time (Schmidt & Patankar 1991; 
Stock & Haase 1999; Walters & Leylek 2004; see also Rurnsey et al. 2005 for a more 
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complete list). These efforts have mainly focused on engineering approaches rather 
than on approaches more theoretically based. The latter approaches have received 
much less attention (e.g., Thacker et al.l999a,b; Jovanovic & Pashtrapanska 2004; 
Rumsey et al.2005) and as yet have not reached a point where extensive flow field 
predictions have been made. It is worth pointing out some characteristics that 
a predictive transition-sensitized turbulence model should possess. These include: 
properly identifying the quantities actually being computed by such models, not 
requiring a priori input of transition or re-laminarization location, and properly 
characterizing the dynamics that such quantities should possess. For methods uti- 
lizing closed transport equations for dynamic variables, such as disturbance kinetic 
energy or disturbance velocity second- moments (stresses), the intent is to compute 
the entire flow domain such as upstream, over and downstream of airfoils, wings, 
turbine blades, or aerodynamic bodies. It is necessary to formulate equations ca- 
pable of describing both viscous dominated and turbulent fluctuations within the 
same flow. 

In general, the usual types of transition modes considered are natural and by- 
pass transition. In natural transition, the fluid responds to a small random dis- 
turbance, which spatially and temporally evolves into a complex nonlinear process 
leading to a fully turbulent flow. The statistical correlations associated with these 
disturbances in the laminar transition region can be described by transport equa- 
tions. While stability theory focuses on the description of deterministic quantities, 
the turbulent prediction methods such as RANS use equations for statistical corre- 
lations. Thus, these transitioning disturbance fields are represented by correlations 
of fluctuating quantities with an associated probability density function (Thacker et 
a!.1999a, 6), and are governed by the same type of correlation transport equations 
as the RANS stresses in the fully turbulent region. In by-pass transition or in a 
re-laminarizing turbulent flow, the flow dynamics are governed by either turbulent 
fluctuations in the free-stream (by-pass transition) or an existing turbulent field be- 
ing subjected to some flow variation (such as pressure-gradients) causing the flow 
to re-laminarize. 

Although the characterizations of the disturbance field just described validate 
the use of transport equations that are similar in form to the RANS equations, 
the commonly used engineering prediction methods have directly applied modified 
RANS models to the prediction of such transitioning flows. The most common and 
simple approach to handling a transition region on an aerodynamic body is to 
control the energy production term by having it “switched off” in the transition 
region. Unfortunately, this a priori knowledge limits the range of applicability and 
general predictive nature of such models. More sophisticated approaches have been 
employed, such as the correlation-based method of Langtry and Menter (2005), 
which is built strictly on local variables and solves two transport equations, one 
for intermittency and one for a transition onset criterion. These currently accepted 
engineering approaches to delimiting the transitioning regions of the flow need to be 
investigated further. Criterion should be developed (e.g. Jovanovic & Pashtrapanska 
2004) that would identify such transitioning zones based on some dynamic measure 
rather than a correlated input invoked by the user. 

Having now identified the type of disturbance correlation fields being computed 
as well as the threshold criterion commonly used delimiting these transitioning 
zones, the next (more subtle) point is whether such correlation fields correctly 


Article submitted to Royal Society 



Modeling Turbulent Aerodynamic Flows 


21 


describe the transitioning or re-laminarization processes. There are two dynamic 
characteristics that the transitional and re-laminarization fields should possess: the 
first is that the “disturbance” or eddy viscosity should be vanishing and second, 
that the ratio of the mean flow time scale (represented by the mean shear S) to 
the disturbance field time scale (represented by the kinetic energy K and kinetic 
energy dissipation rate e) vanish. The former dynamic characteristic is easily met 
and is usually achieved in existing models by “turning off” the production term in 
either the kinetic energy or second-moment equations in the transitioning region. 

Recent studies by Rumsey et a/. (2006) and Petterson-Reif et a/. (2006) have 
identified some distinguishing characteristics of several two-equation linear eddy 
viscosity and algebraic stress models that can have an impact on any attempts 
to utilize such formulations in predicting transitioning and re-laminarizing flows. 
These studies have shown through a dynamical systems analysis using nullclines in 
K and e phase space that solution regions where disturbance (eddy) viscosities given 
by K 2 /e, for example, often do not predict the correct laminar limit in terms of the 
mean flow to turbulent time-scale ratio, which should be e/ SK — > 0 in the laminar 
regime. Such regions showing incorrect behavior were termed “pseudo laminar.” 
These results have a direct implication on the use of such models for the prediction of 
transitioning and re-laminarizing flows. Without proper calibration and evaluation, 
RANS type models can yield solutions qualitatively similar to transitioning flow 
behavior but with the incorrect dynamic characteristics. 

As an example of this anomalous behavior, the flow over an RAE2822 airfoil at 
freestream Mach number M = 0.75, angle-of-attack = 2.72, and Re = 6.2 x 10 6 was 
computed. A plot showing the airfoil shape and resulting pressure contours for these 
conditions is given in Fig. 5. There is a strong shock wave present on the airfoil 
upper surface near 65% chord; whereas, the flow on the lower surface remains sub- 
sonic. In this example, the initial and boundary conditions were kept fixed but the 
numerical solution method was changed (for details see Rumsey et aZ.2006). Figure 
6 shows the skin-friction distribution on the lower surface of the airfoil obtained 
using two different numerical solution strategies for obtaining converged solutions. 
The converged results were completely different, with each suggesting a transition 
location in a different place. This inconsistent behavior was due to the fact that 
particular forms of the K-e model can converge to either a “pseudo-laminar” state 
or the intended turbulent state, depending on initial conditions or other numerical 
parameters. Figure 7 shows the computed results along with the theoretical trajec- 
tories at a point in the boundary layer on a flat plate that eventually converged to 
a “pseudo- laminar” result. Here, diffusion effects were negligible, and there is ex- 
cellent agreement between theory and computation. The SK — e phase plot shows 
that the degenerate fixed point at SK = e = 0 is reached. 

These results indicate that any attempt at using transport equations for statis- 
tical correlations such as kinetic energy and associated dissipation rate can yield 
solutions dependent on the character of the equations rather than on any physical 
calibration. Thus, the utilization of such equations in predicting flows containing 
transitioning and re-laminarizing regions can be problematic. A careful assessment 
of the solutions is needed to determine whether they are replicating the calibration 
characteristics, or are a result of phase plane trajectories inherent in the form of 
the equations. 
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Figure 5. RAE2822 airfoil solution 
showing static pressure contours 



Figure 6. Streamwise variation of skin- 
friction coefficient on RAE2822 airfoil 
lower surface for two different solution 
procedures 



Figure 7. Comparison of theory with computed result at a point in the “pseudo-laminar” 

of the flat plate 
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( b ) Separation and Control 

Another topic of current interest is flow control, particularly for reduction or 
elimination of separation. By using suction, small jets, or actuators on wings, for 
example, it is possible to delay or eliminate the onset of separation. There have been 
recent attempts to assess both the experimental measurement capabilities as well 
as the model prediction capabilities of representative benchmark flows (Rumsey et 
oZ.2006). In such flows, flow field characteristics not well understood or documented 
need to be investigated. These include the effects of fluid injection through synthetic 
jets into quiescent air and cross-flow, and the effects of steady and unsteady fluid 
injection and suction on separated flow regions. A “hump model” configuration 
(Seifert and Pack 2002) exemplifies the latter effect, and is also a realistic configu- 
ration for flow control on wings. Some recent results can be shown highlighting the 
predictive capabilities of current methods. This flow has been studied extensively 
in both steady and time-dependent flow control applications using a wide vari- 
ety of methods that included RANS, hybrid RANS/LES, LES and (underresolved) 
DNS. The cases studied have included no-control, steady suction, and synthetic jet 
(unsteady blowing/suction) control conditions. 

The subsonic flow boundary layer in this case separates at the back side of the 
hump, near 65% chord. With no flow control, the flow reattaches near x/ c = 1.11. 
With either steady suction or synthetic jet flow control applied near the separation 
point, the separation extent can be either reduced or eliminated. Most RANS mod- 
els appear to do a reasonable job predicting the separation point in these cases, but 
all models ranging from one-equation eddy viscosity through full Reynolds stress 
invariably predict too large a separation extent, either in steady-state or in the 
mean. A typical example is shown in Fig. 8, using a K — to Explicit Albebraic 
Stress Model. In this figure, time-averaged mean results are shown for the synthetic 
jet case, with oscillation frequency 138.5 Hz and peak velocity out of the slot of 
26-27 m/s (freestream velocity was 34.6 m/s). In the experiment (Greenblatt et 
o?.2005), the mean flow reattaches near x/c = 1.0; whereas the RANS method 
predicts later reattachment. The reason for the poor RANS behavior is believed 
to be the fact that it dramatically underpredicts the turbulent shear stress in the 
separated region, as shown for typical results in Fig. 9. This underprediction leads 
to too little mixing, and hence a tendency to remain separated too long. 

On the other hand, blended RANS-LES, LES, and DNS appear to do a better 
job predicting the turbulence characteristics in the separated region. By resolving 
the larger eddies (rather than modeling them), much of the turbulence dynamics 
in the separated region can be properly accounted for, yielding earlier reattach- 
ment. However, these methods are generally very time-consuming to compute, and 
under-resolution in time or space can have a significant negative impact on results. 
Nonetheless, recent LES results (Saric et al. 2006) suggest that careful application of 
this method is capable of capturing important effects seen in the hump experiment. 

Other recent studies have focused on circulation control applications using 
Coanda-type jets over flaps (e.g., Lee- Rausch et al. 2006, Shmilovich and Yadlin 
2006) and over circular trailing edges (e.g., Swanson & Rumsey 2006, Slomski et 
oL2006). In both of these types of aerodynamic flows, flow control can greatly en- 
hance the maximum achievable lift. In the case of computing the Coanda flow over 
circular trailing edges, the flowfield particularly regarding where the jet separates 
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x/c 

Figure 8. Long-time-average streamlines for hump model with synthetic jet flow control. 


- has been found to be very sensitive to both numerical parameters as well as to 
turbulence modeling. For example, proper sensitization of the turbulence model to 
streamline curvature effects appears to be very important. Typical examples show- 
ing streamlines around a Coanda surface at the back of an airfoil are shown in 
Fig. 10. The experiment is due to Novak et a/.(Novak et al. 1987). Here, a steady 
jet issues out of a small slot located on the upper surface near x/c = 0.93. Be- 
cause of the Coanda effect, the jet “sticks” to the airfoil surface and draws the flow 
around to just beyond the back of the circular trailing edge, enhancing the circula- 
tion around the body. In this case, a one-equation linear eddy-viscosity turbulence 
model predicts the jet to separate from the Coanda surface too late; whereas, the 
same model with a curvature correction predicts results in good agreement with 
experiment. However, because it is also very difficult to perform high-quality val- 
idation experiments (that provide sufficient information to evaluate and improve 
turbulence models) on this type of configuration, separation flow control remains 
an area of very active research both experimentally and computationally. 


5. Concluding Remarks 

Significant strides have been made in the prediction of turbulent aerodynamic flows 
in the last quarter century. The availability and applicability of commercially avail- 
able numerical solvers, that routinely solve very complex aerodynamic flow fields, 
has led some to believe that the technical discipline has reached full maturity. Un- 
derstanding and ultimately predicting complex turbulent flow fields has and still 
remains as yet an unsolved problem. With each level of success comes increased 
demands on accuracy and solution capabilities that at its limits remains beyond 
current capabilities. 

The topics addressed in this overview are but a brief glimpse into topical areas 
of research that are currently being investigated and which will help extend our 
ability to both understand and predict even more complex aerodynamic flows. 
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Figure 9. Turbulent shear stress profiles at x/c = 0.8 at four times during the oscillation 
cycle, for hump model with synthetic jet flow control. 
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